// ssr.cpp: implementation of the Successive Survivable Routing class.
//
//////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <math.h>

#include "ssr.h"
#include "gaconfig.h"
#include "gatypes.h"
#include "garandom.h"
//#include "BHeapNode.h" // binHeap
//#include "HeapNode.h" // FibHeap


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

ssr::ssr(char * name, int d=0, int c=0, int n=10, 
		 int ld = 1, int fs = 0, int cf=0,int fsc=0,int pr=0,int dc=0, int mt=0,
		 double pb=0.99999)
{
	_name = name; dbg = d; cst = c; nIter = n; 
	load = ld; fseq = fs; cfun = cf; fsce = fsc; pathr =pr; dcase = dc; multi=mt;
	prob = pb;

	const char *DataDir = "data" ;
#ifdef WIN32
	const char *DirSign = "\\" ;
#else
	const char *DirSign = "/" ;
#endif
	const char *InfoFile = "info.dat";	
	const char *TopoFile = "topo.dat";
	const char *HFile    = "H.dat";
	const char *HinfFile = "infoL.dat";
	const char *WorkFile = "work.txt";
	const char *SpaTx = "spare.txt";
	const char *SpaTxND = "spareND.txt";
	const char *SpaTxFD = "spareFD.txt";
	const char *SpaTxFD2 = "spareFD2.txt";
	const char *TrafDat = "traf.dat";
	const char *WorkDat = "work.dat";
	const char *SpaDat = "spare.dat";
	const char *SpaDaND = "spareND.dat";
	const char *SpaDaFD = "spareFD.dat";
	const char *SpaDaFD2 = "spareFD2.dat";
	const char *SPM = "spm.txt";
	const char *SPMND = "spmND.txt";
	const char *SPMFD = "spmFD.txt";
	const char *SPMFD2 = "spmFD2.txt";
	const char *RSV = "rsv.txt";
	const char *RSVND = "rsvND.txt";
	const char *RSVFD = "rsvFD.txt";
	const char *RSVFD2 = "rsvFD2.txt";
	char fullname[200], filename[200];

	strcpy(fullname, "..");
	strcat(fullname, DirSign);
	strcat(fullname, DataDir);
	strcat(fullname, DirSign);
	strcat(fullname, _name);
	strcat(fullname, DirSign);

	ifstream infile;

	// read information of networks
	strcpy(filename, fullname);
	strcat(filename, InfoFile);
	infile.open(filename, ios::in | ios::nocreate);
	if (!infile) {
		printf("Error in finding datafiles...%s\n",filename);
		exit(1);
	}
	infile >> _nNode >> _nLink;
	infile.close();
	if (_nNode > MAX_NODE || _nLink > MAX_LINK){
		printf("This code has limited for %2d nodes and %2d links!\n",
			MAX_NODE, MAX_LINK);
		exit(-1);
	}

	strcpy(filename, fullname);
	strcat(filename, TopoFile);
	infile.open(filename, ios::in | ios::nocreate);
	assert(infile);

	// init adjm and link arrays, adjl
	int i, j;
	for (i=0;i<_nNode;i++){
		adjl[i] = NULL;
		nodelist[i].v = i; nodelist[i].d = MAX_COST;
		for (j=0;j<_nNode;j++)
			adjm[i][j] = -1;
	}
	int n1, n2;
	nodeptr t;
	for (i=0; i < _nLink; i++){
		infile >> n1 >> n2 ;
		n1--; n2--;
		assert(0<=n1 && n1<n2 && n2 <_nNode); //range check
		link[2*i][0] = n1;
		link[2*i][1] = n2;
		link[2*i+1][0] = n2;
		link[2*i+1][1] = n1;
		adjm[n1][n2] = 2*i;
		adjm[n2][n1] = 2*i + 1;
		t = (nodeptr)malloc(sizeof(struct node_str));
		t->v = n2; t->next = adjl[n1]; adjl[n1] = t;
		
		t = (nodeptr)malloc(sizeof(struct node_str));
		t->v = n1; t->next = adjl[n2]; adjl[n2] = t; 
	}
	infile.close();
	//if (dbg) 
	//	cout << "Topology input finished.\n";

	memset(F, 0, sizeof(char)*MAX_FAIL*MAX_LINK*2);
	memset(LR, 0, sizeof(char)*MAX_LINK*MAX_LINK*4);
	memset(U, 0, sizeof(char)*MAX_DMD*MAX_LINK*4); //U[MAX_DMD*2][MAX_LINK*2]
	memset(Atabu, 0, sizeof(char)*MAX_DMD*MAX_LINK*4); //B[2][MAX_DMD*2][MAX_LINK*2]
	switch (fsce) {
	case 0: 
		_nFail = _nLink;
		for (i=0; i < _nLink; i++){
			unitcost[2*i] = 1;
			unitcost[2*i+1] = 1;
			F[i][2*i] = 1; 
			F[i][2*i+1] = 1;
		}
		break;
	case 1:
		_nFail = _nNode;
		for (i=0; i < _nLink; i++){
			unitcost[2*i] = 1;
			unitcost[2*i+1] = 1;
			n1 = link[2*i][0];
			n2 = link[2*i][1];
			F[n1][2*i] = 1; F[n1][i*2+1] = 1;
			F[n2][2*i] = 1; F[n2][i*2+1] = 1;
		}
		break;
	default:
		cout << " Not implemented to read failure scenarios..."<<endl;
		assert(1);
	}
	if (multi > 0) {
		strcpy(filename, fullname);
		strcat(filename, HinfFile);
		infile.open(filename, ios::in | ios::nocreate);
		assert(infile);
		infile >> n1 >> n2;
		_nlNode = n1;
		_nlLink = n2;
		infile.close();

		strcpy(filename, fullname);
		strcat(filename, HFile);
		infile.open(filename, ios::in | ios::nocreate);
		assert(infile);
		for (i=0; i < _nLink; i++){
			//if (multi!=3){
				unitcost[2*i] = 0;
				unitcost[2*i+1] = 0;
			//}
			for (j=0; j<_nlLink;j++){  // F = H^T
				infile >> n1;
				H[2*i][j] = n1;
				H[2*i+1][j] = n1;
				//if (multi!=3){
					unitcost[2*i] += n1;
					unitcost[2*i+1] += n1;
				//}
			}
		}
		infile.close();
	}
	if  (multi >=2) { // consider H in optimization
		_nFail = _nlLink;  // number of links at lower layer
		for (i=0; i < _nLink; i++){
			for (j=0; j<_nFail;j++){  // F = H^T
				F[j][2*i] = H[2*i][j];
				F[j][2*i+1] = H[2*i+1][j];
			}
		}
	}

	// init flows (demand pairs)
	strcpy(filename, fullname);
	strcat(filename, TrafDat);
	outfile6.open(filename, ios::out);
	assert(outfile6);

	_nDmd = 0; // use full meshed load or all more-than-one-hop load.
	for (i = 0; i<_nNode-1; i++){
		dmdm[i][i] = -1;
		for (j=i+1; j<_nNode; j++) {
			if (dcase ==0 || (dcase==1 && adjm[i][j]!= -1)){
				flow[_nDmd][0] = i;	flow[_nDmd][1] = j;
				dmdm[i][j] = _nDmd; _nDmd ++;
				flow[_nDmd][0] = j;	flow[_nDmd][1] = i;
				dmdm[j][i] = _nDmd; _nDmd ++;
			}else {
				dmdm[i][j] = -1; dmdm[j][i] = -1;
			}
		}
	}
	dmdm[_nNode-1][_nNode-1] = -1;
	
	for (i = 0; i<_nNode; i++){
		for (j=0; j<_nNode; j++) {
			if (i == j)
				outfile6 << " ";
			else
				outfile6 << M(dmdm[i][j]) <<" ";
		}
		outfile6 <<endl;
	}
	
	strcpy(filename, fullname);
	strcat(filename, WorkFile);
	outfile1.open(filename, ios::out);
	//outfile1.setmode(filebuf::text);
	assert(outfile1);
	strcpy(filename, fullname);
	if (pathr == 1)
		strcat(filename, SpaTxFD);
	else if (pathr == 2)
		strcat(filename, SpaTxFD2);
	else if (pathr == 3)
		strcat(filename, SpaTxFD);
	else if (fsce == 1)
		strcat(filename, SpaTxND);
	else
		strcat(filename, SpaTx);
	outfile2.open(filename, ios::out);
	//outfile2.setmode(filebuf::text);
	assert(outfile2);
	strcpy(filename, fullname);
	if (pathr == 1 || pathr == 3)
		strcat(filename, SPMFD);
	else if (pathr == 2 || pathr == 4)
		strcat(filename, SPMFD2);
	else if (fsce == 1)
		strcat(filename, SPMND);
	else
		strcat(filename, SPM);
	outfile3.open(filename, ios::out);
	//outfile3.setmode(filebuf::text);
	assert(outfile3);
	strcpy(filename, fullname);
	if (pathr == 1 || pathr == 3)
		strcat(filename, RSVFD);
	else if (pathr == 2|| pathr == 4)
		strcat(filename, RSVFD2);
	else if (fsce == 1)
		strcat(filename, RSVND);
	else
		strcat(filename, RSV);
	outfile4.open(filename, ios::out);
	assert(outfile4);
	strcpy(filename, fullname);
	strcat(filename, WorkDat);
	outfile5.open(filename, ios::out);
	assert(outfile5);

	strcpy(filename, fullname);
	if (pathr == 1 || pathr == 3)
		strcat(filename, SpaDaFD);
	else if (pathr == 2|| pathr == 4)
		strcat(filename, SpaDaFD2);
	else if (fsce == 1)
		strcat(filename, SpaDaND);
	else
		strcat(filename, SpaDat);
	outfile51.open(filename, ios::out);
	assert(outfile51);

	workBW = 0;	
	memset(A, 0, sizeof(char)*MAX_DMD*MAX_LINK*4); //A[MAX_DMD*2][MAX_LINK*2]
	memset(hopA, 0, sizeof(int)*MAX_DMD*2);
	memset(MF,0, sizeof(int)*MAX_DMD);
	memset(Aroute, -1 ,sizeof(char)*MAX_DMD*2*MAX_NODE); 
	memset(Artlink, -1 ,sizeof(char)*2*MAX_DMD*MAX_NODE);
	memset(Artnode, -1 ,sizeof(char)*2*MAX_DMD*MAX_NODE);
	
	// dangerous line here!!!
	// when number of links is larger than 128

	NLC[0] = 0;NLC[1] = 3.7;
	for (i=2;i<2000;i++) {
		if (i <=3)
			NLC[i] = 7.72; //3.7 * pow(3,0.67);
		else if (i <=12)
			NLC[i] = 19.55; //3.7 * pow(12,0.67);
		else if (i <=48)
			NLC[i] = 49.50; //3.7 * pow(48,0.67);
		else if (i <=192)
			NLC[i] = 125.32; //3.7 * pow(192,0.67);
		else
			NLC[i] = (1+i/192) * NLC[192];
	}

	// initialize the link availability 
	// since we used single value for all links, 
	// there is no need for link-wise values.

	init();
}

ssr::~ssr() 
{
	outfile1.close();	outfile2.close();
	outfile3.close();	outfile4.close();
	outfile5.close();	outfile51.close();
	outfile6.close();
}

int ssr::init()
{
	// clear backup paths
#ifndef FAILDENP
	memset(B, 0, sizeof(char)*MAX_DMD*MAX_LINK*4); 
	memset(hop, 0, sizeof(int)*MAX_DMD*2);
	memset(route, -1 ,sizeof(char)*2*MAX_DMD*MAX_NODE);
	memset(rtlink, -1 ,sizeof(char)*2*MAX_DMD*MAX_NODE);
#else
	memset(B, 0, sizeof(char)*MAX_DMD*MAX_LINK*MAX_LINK*8);
	memset(hop, 0, sizeof(int)*MAX_LINK*MAX_DMD*4);
	memset(route, -1 ,sizeof(char)*MAX_LINK*4*MAX_DMD*MAX_NODE);
	memset(rtlink, -1 ,sizeof(char)*MAX_LINK*4*MAX_DMD*MAX_NODE);
#endif

	spareBW = 0;
	memset(C, 0, sizeof(int)*MAX_LINK*MAX_LINK*4);
	memset(C2, 0, sizeof(int)*MAX_LINK*MAX_LINK);
	memset(rsvBW, 0, sizeof(int)*MAX_LINK*6);
	return 0;
}

int ssr::visit(int k)
{
	static int id;
	nodeptr t;
	int m, mintop;  // these two variable are added to deal with 2-connectivity

	id ++;
	val[k] = id; mintop = id+1;
	t = adjl[k];
	while (t != NULL) {
		if (val[t->v] == 0) {
			use[adjm[k][t->v]] = 1; // mark this link has been used
			m = visit(t->v);
			if (m < mintop) 
				mintop = m;
			if (m > val[k]) {
				conn2 = (adjm[k][t->v] % _nLink) +1;
				printf("Warning: criticle link %2d between nodes %2d and %2d !\n",
					conn2, t->v +1, k+1);
			}
		}
		else if (val[t->v] <mintop && use[adjm[t->v][k]] ==0 )
			mintop = val[t->v];

		t = t->next;
	}
	return mintop;
}

int ssr::listdfs(void)
{
	memset(val,0,sizeof(int)*MAX_NODE);
	memset(use,0,sizeof(int)*MAX_LINK*2);
	conn2 = 0; // when it is turned to none zero, broke some link will seperate the graph
	int conn = visit(0);  // return at least 1.  larger than 1 means not 2-link-connected
	for (int k=0; k<_nNode; k++){
		if (val[k] ==0) {
			printf("Warning: node %2d is not reachable from node 1 !\n",k+1);
			conn = 0;
		}
	}

	if (conn == 0)
		printf("The graph is not connected!\n"); 
	else {
		printf("Network %s is connected ", _name);
		if (conn2 == 0 && conn == 1) {
			printf("and two-link-connected.\n");
			conn = 2;
		}
		else {
			printf("BUT not 2-link-connected.\n");
			conn = 1;
		}
	}	
	return conn;
}

// should change to route instead of B here.....!!!!!
int ssr::addC(int dmd)
{
	int temp;
	for (int k=0; k<_nFail; k++){
		if (U[dmd][k] == 1 ) {
			for (int m=0;m<_nLink; m++){
				if (B[0][dmd][m*2]==1 || B[0][dmd][m*2+1]==1) {
					C[m*2][k] += M(dmd); C[m*2+1][k] += M(dmd);
	}}}}
	if (multi==3){
		for (int k=0; k<_nFail; k++){
			if (U[dmd][k] == 1 ) {
				for (int m=0;m<_nlLink; m++){
					temp = 0;
					for (int j=0;j<_nLink*2; j++)
						temp += H[j][m]*B[0][dmd][j];
					C2[m][k] += (temp*M(dmd)*2);
				}
			}
		}
	}

	getrsvBW();
	return 0;
}

int ssr::minusC(int dmd)
{
	int temp;
	for (int k=0; k<_nFail; k++){
		if (U[dmd][k] ==1) { 
			for (int m=0;m<_nLink; m++){
				if (B[0][dmd][m*2]==1 || B[0][dmd][m*2+1]==1 ) {
					C[m*2][k] -= M(dmd); C[m*2+1][k] -= M(dmd);
	}}}}
	if (multi==3){
		for (int k=0; k<_nFail; k++){
			if (U[dmd][k] == 1 ) {
				for (int m=0;m<_nlLink; m++){
					temp = 0;
					for (int j=0;j<_nLink*2; j++)
						temp += H[j][m]*B[0][dmd][j];
					C2[m][k] -= (temp*M(dmd)*2);
				}
			}
		}
	}
	getrsvBW();
	return 0;
}


// given C, find cost on link for demand 

void ssr::getrsvBW(void)
{
	int k,j,m, temp;
	for (m=0; m<_nLink*2; m++){
		temp = 0;
		for (k = 0; k <_nFail; k++)
			temp = max(temp, C[m][k]);
		rsvBW[1][m] = temp; 
	}
	if (multi==3) { //s = H^T * C  or G
		for (j=0; j<_nlLink; j++){
			temp = 0; 
			for (k=0; k<_nFail; k++)
				temp = max(temp, C2[j][k]);
			rsvBW[2][j] = temp; 
		}
	}
}

// given src/dst and cost, return min cost and path in succ.
// need: link, cost, _nNode, _nLink;
int ssr::spf1(int dmd)
{  // looks like a BFM algorithm ??

	int k,i,j;
	double dtemp;
	double d[MAX_NODE];
	int dnew[MAX_NODE];
	
	int src = flow[dmd][0];
	int dst = flow[dmd][1];
	memset(dnew,0,sizeof(int)*MAX_NODE);
	for (i=0; i < _nNode; i++){
		d[i]=MAX_COST;
		succ[i]=-1;
	}
	/*for (i = 0; i<_nNode; i++){
		d[i] = MAX_NODE;
		succ[i] = -1;
	}*/
	d[dst] = 0;
	dnew[dst] = 1;

	do{
		for (k = 0; k<_nLink*2; k++){
			i=link[k][0]; j=link[k][1];
			if (d[i] > d[j] + cost[k]) {
				d[i] = d[j] + cost[k];
				succ[i] = j;
			}
		}

		// Find a node with the smallest label i
		dtemp = MAX_COST;
		i = MAX_NODE;
		for (j=0; j <_nNode; j++){
			if (dnew[j] ==0 && dtemp > d[j]) {
				dtemp = d[j]; //used for hop count here
				i = j;
			}
		}
		if (i != MAX_NODE)
			dnew[i] = 1; // fix this node i
		else 
			return(-1); // if no update, it means there is no solution
	}while (dnew[src] != 1);
	return 0;  // this is not the min hop anymore.  It is min cost
}


// return hop count of the minimum path.
int ssr::findhop(int dmd)
{	
	int k,j,i = flow[dmd][0];
	int dst = flow[dmd][1];
	int dtemp = 0;
	memset(path,0,sizeof(char)*MAX_LINK*2);
	do {
		j = succ[i];
		k = adjm[i][j];
		path[k] = 1;
		i=j;
		dtemp ++;
	}while ( j != dst);
	//assert(i==dst); // impossible state
	return dtemp;  // this is the hop count
}

// augmenting path algorithm
int ssr::apa(int dmd)
{
	int k2, hp;
	//try to find the second link disjoint path.

	do {
		int k;
		findworking(dmd);

		for (k = 0; k < _nLink*2 ; k ++){
			if (Atabu[dmd][k] == 1 || A[dmd][k]==-1) 
				cost[k] = MAX_COST;
			else
				cost[k] = 1;
		}

		// check if current working path form a trap topology, using augmenting path algo.
		if ( spf(dmd) != -1){
			for (int k = 0; k < _nLink*2 ; k ++){
				if (A[dmd][k] == -1)
					A[dmd][k] = 0;
			}
			// clean trap links.
			return 0; // find a backup path, exit apa.
		}
		// reverse direction of working path is enabled 
		// only disjointed on directed links.
		for (k=0; k<_nLink*2; k+=2){
			if (A[dmd][k] == 1)
				cost[k+1] = -0.99;
			else if (A[dmd][k+1] == 1)
				cost[k] = -0.99;  // reverse is enabled
			else if (A[dmd][k+1]!=-1) {
				cost[k] = 1; cost[k+1]=1;
			}
		}

		if (fsce ==1 ) 
			expandnode(dmd);
		hp = spf(dmd);
		if ( hp  != -1){
			if ( fsce ==1 ) 
				contractnode(dmd);	

			hopA[dmd] = findhop(dmd);  // get path[] from succ[].
			//  now we have second path in augmenting path algorithm
			// find trap link
			for (k=0; k<_nLink*2; k+=2) {
				k2 = k+1;
				if (A[dmd][k]==1 && path[k2] ==1) {
					A[dmd][k] = -1; A[dmd][k2] = -1;// marked as trap link
					cout << "trap for" << flow[dmd][0] <<"-"<<flow[dmd][1];
					cout << ": " << link[k][0] <<"-"<<link[k][1]<<endl;
				} else if (A[dmd][k2]==1 && path[k] ==1) {
					A[dmd][k2] = -1; A[dmd][k] = -1;// marked as trap link
						cout << "trap for" << flow[dmd][0] <<"-"<<flow[dmd][1];
					cout << ": " << link[k2][0] <<"-"<<link[k2][1]<<endl;
				}
				if (A[dmd][k] !=-1) 
					A[dmd][k] = 0; 
				if (A[dmd][k2] !=-1) 
					A[dmd][k2] = 0;  // keep trap link for next turn
				/*else if (A[dmd][k2]==1 || path[k] ==1 ||
				     A[dmd][k]==1 || path[k2] ==1){
				A[dmd][k] = 0; A[dmd][k2] = 0;
				} else {
					A[dmd][k] = -1; A[dmd][k2] = -1;
				}
				//this part is a short cut but can not get min working.
				*/	
			}  //for 	
			hopA[dmd] = 0;  // prepare for new working path
			//find working path again.	
		} else {
			printf("Wrong position!\n"); exit(1);
		}	
	}while (1);
}

int ssr::expandnode(int dmd)
{ // this process build a residual topology after split path nodes.
	int i,j,j2;
	nodeptr t;
	int clink = _nLink*2, cnode = _nNode; 
	slink = _nLink; snode = _nNode; 
	
	for (i = 0; i<_nFail; i++){
		if (U[dmd][i] ==1){ 
			// split this node, add one more link
			link[clink][1] = i;
			link[clink][0] = cnode;
			cost[clink] = 0; // allow routing
			clink++; 
			link[clink][0] = i;
			link[clink][1] = cnode;
			cost[clink] = MAX_COST; // not allow routing
			clink++; 
			smap[cnode] = i;
			cnode ++;
			// for all outgoing li
			for (t=adjl[i]; t!=NULL;t=t->next){
				j = adjm[i][t->v]; j2=adjm[t->v][i]; 
				// link from n1 to t->v is adjm[n1][t->v] 
				// move outbound link to splited node when 
				// it is not on the path
				if (A[dmd][j] !=1 && A[dmd][j2] !=1)
					link[j][0] = cnode -1;
				if (A[dmd][j] == 1)
					link[j2][1] = cnode -1;
			}
		}
	}
	_nLink = clink/2; _nNode = cnode;
	return 0;
}
int ssr::contractnode(int dmd)
{
	int i,j,j2;
	nodeptr t;
	int clink = _nLink*2, cnode = _nNode; 
	_nLink = slink; _nNode = snode ; 
	
	for (i=flow[dmd][0]; succ[i]!=flow[dmd][1]; i=succ[i]){
		if (succ[i] >= _nNode){
			if (succ[succ[i]] != smap[succ[i]])
				succ[smap[succ[i]]] = succ[succ[i]];
			succ[i] = smap[succ[i]]; // jump through vitual node
		}
	}

	for (i = 0; i<_nFail; i++){
		if (U[dmd][i] ==1){ 
			// put all outgoing link back to i
			for (t=adjl[i]; t!=NULL;t=t->next){
				j = adjm[i][t->v]; 
				j2=adjm[t->v][i]; 
				if (A[dmd][j] !=1 && A[dmd][j2] !=1)
					link[j][0] = i;
				if (A[dmd][j] == 1)
					link[j2][1] = i;
			}
		}
	}
	return 0;
}
//inline 
double ssr::linkCost(int link, int dmd)
{
	if (cst == 1 || cst ==3) return 1; // For RAFT and NS
	double augRsv = 0;
	int t;
	// using spm matrix C
	// link %= _nLink;
	int current = rsvBW[1][link];
	int curTotal = rsvBW[0][link]+rsvBW[1][link];
	int augtemp ;
	if (multi==3) { //s = H^T * C  , v^l first, then v = H*v^l
		for (int l=0;l<_nlLink;l++) {
			augtemp =0;	
			for (int t=0; t<_nFail; t++){ //find v^l
				if (U[dmd][t]==1) {
					augtemp = max(augtemp,C2[l][t] + M(dmd) - rsvBW[2][l]);
				}
			}
			if (H[link][l] == 1) 
				augRsv += augtemp;
		}
		return augRsv;
	}
	if (cst == 0){
		for (t=0; t<_nFail; t++){
			if (U[dmd][t]==1)
				augRsv = max(augRsv, C[link][t] + M(dmd) - current );
		}
	}
	else {//cst ==2;  // Kodialam's cost function.
		augRsv = max(0,MF[dmd] + M(dmd) - current);
	}
	// if unit cost of link is used
	augRsv *= unitcost[link];  
	if (cfun ==1 ) //if non-linear, return additional cost
		augRsv = NLC[(int)ceil(augRsv)+curTotal] - NLC[curTotal];
		
	return augRsv;
}

int ssr::findbackup(int dmd)
{
	int src = flow[dmd][0];
	int dst = flow[dmd][1];
	//int dmd2 = dmdm[dst][src];
	int i,k2;

	if (hop[0][dmd]>0){  // only if this flow has backup
		//removebackup(dmd); removebackup(dmd2);
		minusC(dmd); 
		//minusC(dmd2);
	}
	
	for (int k = 0; k < _nLink*2 ; k ++){
		if (Atabu[dmd][k] == 1) 
			cost[k] = MAX_COST;
		else
			cost[k] = linkCost(k, dmd);
	}

	k2 = spf(dmd);
	if ( k2 != -1){ // find SP
		hop[0][dmd] = findhop(dmd);	
 
		if (memcmp(path, B[0][dmd],sizeof(char)*MAX_LINK*2)){
		// is updated
			memcpy(B[0][dmd],path,sizeof(char)*MAX_LINK*2);
			memcpy(route[0][dmd],succ,sizeof(char)*MAX_NODE);
			succ2rt(dmd, rtnode[0][dmd], succ);
			Ccurrent ++; 
			k2 = 1;

			for (i=0; i<_nLink;i++){
				B[0][dmd+1][i*2] = B[0][dmd][i*2+1];
				B[0][dmd+1][i*2+1] = B[0][dmd][i*2];
			}	
			hop[0][dmd+1] = hop[0][dmd];
			for (i=0; i<=hop[0][dmd]; i++)
				rtnode[0][dmd+1][i] = rtnode[0][dmd][hop[0][dmd]-i];

			for (i=0; i<_nNode; i++) 
				route[0][dmd+1][i] = -1;
			for (i=0; i<hop[0][dmd]; i++)
				route[0][dmd+1][rtnode[0][dmd+1][i]] 
				= rtnode[0][dmd+1][i+1];


		} else	
			k2 = 0;
  
		Cflag[dmd] = Ccurrent;
		addC(dmd); 
	} else {
		if (fsce ==0) {
			printf("This state is wrong...\n"); 
			exit(-1); 
		} else {
			memset(B[0][dmd],0,sizeof(char)*MAX_LINK*2);
			k2 = 0;
		}
	}
	return k2;

}

#ifdef FAILDENP
int ssr::addCFD(int dmd, int k)
{
	//k = k ;
	if (U[dmd][k] != 1 ) 
		assert(1);
	for (int i = dmd; i <dmd+1; i++){
		for (int m=0; m<_nLink; m++){
			if (B[k][i][m*2]==1 || B[k][i][m*2+1]==1) {
				C[m*2][k] += M(i); C[m*2+1][k] += M(i);
			}
			if (B[k][i][m*2]==-1 || B[k][i][m*2+1]==-1){
				C[m*2][k] -= M(i); C[m*2+1][k] -= M(i);
			}
		}
	}
	getrsvBW();
	return 0;
}

int ssr::minusCFD(int dmd, int k)
{
	if (U[dmd][k] != 1 ) 
		assert(0);
	for (int i = dmd; i <dmd+1; i++){
		for (int m=0;m<_nLink; m++){
			if (B[k][i][m*2]==1 || B[k][i][m*2+1]==1) {
				C[m*2][k] -= M(i); C[m*2+1][k] -= M(i);
			}
			if (B[k][i][m*2]==-1 || B[k][i][m*2+1]==-1) {
				C[m*2][k] += M(i);C[m*2+1][k] += M(i);
			}
		}
	}
	getrsvBW();
	return 0;
}

double ssr::linkCostFD(int link, int dmd, int j)
{
	if (cst == 1 || cst ==3) return 1; // For RAFT and NS
	double augRsv = 0;
	int t;
	// using spm matrix C
	//link %= _nLink;
	int current = rsvBW[1][link];
	int curTotal = rsvBW[0][link]+rsvBW[1][link];
	if (cst == 0){
		t=j;
		if (U[dmd][t]==1)
			augRsv = max(augRsv, C[link][t] +M(dmd) - current );
	}
	else {//cst ==2;  // Kodialam's cost function.
		augRsv = max(augRsv,MF[dmd] + M(dmd) - current);
	}
	// if unit cost of link is used
	augRsv *= unitcost[link];  
	if (cfun ==1 ) //if non-linear, return additional cost
		augRsv = NLC[(int)ceil(augRsv)+curTotal] - NLC[curTotal];
		
	return augRsv;
}

int ssr::stub_release(int dmd, int j)
{
	int i, s, t, k;

	for (i=0; i< hopA[dmd]; i++){
		s = Artnode[dmd][i];
		t = Artnode[dmd][i+1];
		k = adjm[s][t];
		if (k < 0) 
			assert(0); // error in working route.
		if (F[j][k] != 1 )
			path[k] --;
	}
	return 0;
}

int ssr::findbackupFD(int dmd)
{
	int src = flow[dmd][0];
	int dst = flow[dmd][1];
	int i, j, k, k2;

	for (j = 0; j <_nFail; j++) {
	if (U[dmd][j]==1){

		if (hop[j][dmd]>0)
			minusCFD(dmd, j); 

		for (k= 0; k < _nLink ; k ++){
			if (A[dmd][k*2] == 1 || A[dmd][k*2+1] == 1) {
				cost[k*2] = MAX_COST;cost[k*2+1] = cost[k*2] ;
			} else {
				cost[k*2] = linkCostFD(k*2, dmd, j);
				cost[k*2+1] = linkCostFD(k*2+1, dmd, j);
			}
		}
		k2 = spf(dmd); 
		if ( k2 != -1){ // find SP
			hop[j][dmd] = findhop(dmd);	 //get path from succ
			if (pathr == 2 || pathr==4) 
				stub_release(dmd, j); // for link j failure
			
			if (memcmp(path, B[j][dmd],sizeof(char)*MAX_LINK*2)){
				// is updated
				memcpy(B[j][dmd],path,sizeof(char)*MAX_LINK*2);
				memcpy(route[j][dmd],succ,sizeof(char)*MAX_NODE);
				succ2rt(dmd, rtnode[j][dmd], succ);
				Ccurrent ++; 
				k2 = 1;

				for (i=0; i<_nLink;i++){
					B[j][dmd+1][i*2] = B[j][dmd][i*2+1];
					B[j][dmd+1][i*2+1] = B[j][dmd][i*2];
				}	

				hop[j][dmd+1] = hop[j][dmd];
				for (i=0; i<=hop[j][dmd]; i++)
					rtnode[j][dmd+1][i] = rtnode[j][dmd][hop[j][dmd]-i];

				for (i=0; i<_nNode; i++) 
					route[j][dmd+1][i] = -1;
				for (i=0; i<hop[j][dmd]; i++)
					route[j][dmd+1][rtnode[j][dmd+1][i]] 
					= rtnode[j][dmd+1][i+1];

				Cflag[dmd] = Ccurrent;
			} else 
				k2 = 0;
			addCFD(dmd, j); 
		} else {
				printf("This state is wrong...\n"); 
				exit(-1); 
		}
	} // if (U==1
	} // for (j
	return k2;
}

#endif

// make symetric paths link reversed
int ssr::symetry(int t, int dmd)
{
	assert(t==0 || t==1);
	int dmd2=dmd+ _nDmd, k2;
	//memcpy(B[t][dmd],B[t][dmd2]*2+1,sizeof(char)*_nLink);
	//memcpy(B[t][dmd]*2+1,B[t][dmd2],sizeof(char)*_nLink);
	for (int k =0; k <_nLink; k++){
		k2 = k*2+1;
		B[t][dmd2][k2] = B[t][dmd][k];
		B[t][dmd2][k] = B[t][dmd][k2];
	}
	hop[t][dmd2]= hop[t][dmd];

	int i=flow[dmd][0], j;
	int dst=flow[dmd][1];
	do {
		j=route[t][dmd][i];
		route[t][dmd2][j] = i;
		i=j;
	} while (j != dst);
	
	return 0;
}


// Find Working path for demand and its symetric demand
int ssr::findworking(int dmd)
{
	int k,i,j;
	int src = flow[dmd][0];
	int dst = flow[dmd][1];
	//int dmd2 = dmdm[dst][src];

	for (k = 0; k < _nLink*2 ; k ++){
		cost[k] = 1;
		if (A[dmd][k] == -1){  // used for AugPathAlgo to avoid trap link
			cost[k] = MAX_COST;
			// A[dmd][k] = 0;
		}
	}

	// given src/dst,link, and cost, return int[MAX_LINK] as the shortest path.
	spf(dmd);

	hopA[dmd] = findhop(dmd);
	//memcpy(A[dmd],path,sizeof(char)*MAX_LINK*2);
	for (k = 0; k < _nLink*2 ; k ++){
		if (A[dmd][k] != -1)
			A[dmd][k] = path[k];
	}  // keep old -1 to get rid of all trap links
	memcpy(Aroute[dmd],succ,sizeof(char)*MAX_NODE);
	succ2rt(dmd, Artnode[dmd], succ);
	
	// symetry copy
	for (i=0; i<_nLink;i++){
		A[dmd+1][i*2] = A[dmd][i*2+1];
		A[dmd+1][i*2+1] = A[dmd][i*2];
	} 
	hopA[dmd+1] = hopA[dmd];
	for (i=0; i<=hopA[dmd]; i++)
		Artnode[dmd+1][i] = Artnode[dmd][hopA[dmd]-i];
	
	for (i=0; i<_nNode; i++) 
		Aroute[dmd+1][i] = -1;
	for (i=0; i<hopA[dmd]; i++)
		Aroute[dmd+1][Artnode[dmd+1][i]] = Artnode[dmd+1][i+1];
	

	// calculate U[dmd]
	for (i=dmd; i < dmd+2; i++) {
		for (k=0; k<_nFail; k++) {
			U[i][k] = 0;
			// when it is node failure, remove o/d nodes.
			// equiv. U = A \binarytimes F^T - D_o - D_d 
			if (fsce != 1 || (flow[i][0]!=k && flow[i][1]!=k)){
				//U = A \binarytimes F^T
				for (j = 0; j<_nLink*2; j++){
					if (A[i][j] ==1 && F[k][j]==1){
						U[i][k] = 1;
						break;// for \binarytimes instead of \times
					}	
				}
			}
		}
		// find tabu vector for all flows
		for (int j=0; j<_nLink*2;j++){
			Atabu[i][j] = 0; 
			if ( A[i][j]==1)
				Atabu[i][j] = 1; 
			// Atabu = U \binarytimes F
			for (k = 0; k<_nFail;k++){
				if (U[i][k] == 1 && F[k][j] ==1){
					Atabu[i][j] = 1; 
					break;
				}
			}
		}
	}

	return 0;
}

void ssr::succ2rt(int dmd, char *rtnd, const char *succ)
{
	int i = flow[dmd][0], j = 0; 
	int dst = flow[dmd][1];
	rtnd[j] = i; //src
	if (dbg) 
		cout << "Demand " << i << " -> " << dst << ": Working:" << i ;
	do {
		i=succ[i]; j ++;
		rtnd[j] = i;
		//routelink[dmd][j-1] = adjm[route[dmd][j-1]][i];
		//indDmdLnk[dmd][routelink[dmd][j-1]][0] = 1;
		//rsvBW[routelink[dmd][j-1]][0] ++; workBW ++;
		if (dbg) 
			cout << "-" << i;
	}while (i != dst);
}

int ssr::allworking(void)
{
	int k, m;
	// sepecial treat for the trap topology,  
	// can not removed if do not modify working paths
	for (k = 0; k<_nDmd; k+=2) 
		apa(k);
	
	workBW = 0;
	for (k=0; k<_nDmd; k++)
		workBW += ( hopA[k]*M(k));
	
	int temp, temp2=0;
	for (m = 0; m <_nLink*2;m++) {
		temp = 0;
		for (int k=0; k<_nDmd; k++){
			if ( A[k][m] == 1 )
				temp += M(k);
		}
		rsvBW[0][m] = temp;
		temp2 += (temp*unitcost[m]);
	}
	//assert(workBW == temp2);
	workBW = temp2;
	//MF is for Kodialem's SPI use
	for (k=0; k<_nDmd; k++){
		for (m = 0; m<_nLink;m++){
			 if ((A[k][m*2] == 1||A[k][m*2+1] == 1)
				 && (MF[k] < rsvBW[0][m*2]))  
				 MF[k] = rsvBW[0][m*2];
		}
	}

	return workBW ;
}

void ssr::writeworkdat(void)
{
	int i, j, k, m, n;
	for (i = 0; i<_nNode; i++)
	{
		for (j=0; j<_nNode; j++) {
			if (i!=j){
				k = dmdm[j][i];
				n = hopA[k];
				//if (i > j)
					for (m=0; m<=n;m++)
						outfile5 << Artnode[k][m]+1 << " ";
				//else
				//	for (m=0; m<=n;m++)
				//		outfile5 << Artnode[k][n-m]+1 << " ";
				outfile5 << endl;
			}
		}
		outfile5 << endl;
	}
}

void ssr::writesparedat(void)
{
	int i, j, k, m, n, p;
	for (i = 0; i<_nNode-1; i++){
		for (j=i+1; j<_nNode; j++) {
			k = dmdm[i][j];
			if (pathr == 0){
				n = hop[0][k];
				outfile51 << " ";
				for (m=0; m<=hopA[k];m++)
					outfile51 << Artnode[k][m]+1 << " ";
				outfile51 << "\n";
				for (m=0; m<=n;m++)
					outfile51 << rtnode[0][k][m]+1 << " ";
				outfile51 << endl;
			} else { // for path restoration
				outfile51 << " ";
				for (m=0; m<=hopA[k];m++)
					outfile51 << Artnode[k][m]+1 << " ";
				outfile51 << "\n";
				for (p=0; p< _nFail; p++){
					if (U[k][p]==1){
						if (fsce == 0)
							outfile51 << "link "<< link[p*2][0]+1
								<<"-"<<link[p*2][1]+1<<": ";
						else if (fsce==1)
							outfile51 << "node "<<p<<": ";
						n = hop[p][k];
						for (m=0; m<=n;m++)
							outfile51 << rtnode[p][k][m]+1 << " ";
						outfile51 << endl;
					}
				}
				outfile51 << endl;
			}

		}
		//outfile51 << endl;
	}
}

void ssr::writealltxt(void)
{
	int k, m, j;
	for (k = 0; k<_nDmd; k+=2) {
		for (m=0; m<_nLink; m++)
			outfile1 << (int)(A[k][m*2]+A[k][m*2+1])<< " " ; // min node is 1 not 0
		outfile1 << endl;

		if (pathr == 0){
			for (m=0; m<_nLink; m++)
				outfile2 << (int)(B[0][k][m*2]+B[0][k][m*2+1])<< " " ; // min node is 1 not 0
			outfile2 << endl;
		} 
		else { //pathr == 1;
			outfile2 << k+1 << endl;
			for (j = 0; j<_nLink; j++){
				for (m=0; m<_nFail; m++)
					outfile2 << (int)(B[m][k][j*2]+B[m][k][j*2+1])<< " " ;
				outfile2 << endl;
			}
			//outfile2 << endl;
		}
	}
	for (k = 0; k<_nLink; k++) {
		for (m=0; m<_nFail; m++)
			outfile3 << (C[2*k][m])<<" ";
		//+C[k][m*2+1]+C[k*2+1][m]+C[k*2+1][m*2+1]) <<" " ;
		outfile3 << endl;

		outfile4 << rsvBW[0][2*k] << " " << rsvBW[1][2*k] <<endl;
		//assert(rsvBW[1][k] == rsvBW[1][k*2+1]);
	}
	if (multi ==3 ) {
		for (m=0; m<_nFail; m++)
			outfile4 << rsvBW[2][m] << endl;
	}	
}

int ssr::reshuffle(void)
{
	int n1, n2, n3;
	for (int k = 0; k< (5*_nDmd/2); k++) {// shuffle 5 times.  Assume it is enough for randomness
		n1 = GARandomInt(0,_nDmd/2-1)*2;
		n2 = GARandomInt(0,_nDmd/2-1)*2;
		n3 = seq[n1];
		seq[n1]=seq[n2];
		seq[n2]=n3;
	}
	return 0;
}

int ssr::resort(int fseq)
{
	int n1;
	// sort according to hop count
	if (fseq == 1 || fseq == -1 ||
		fseq == 3 || fseq == -3||
		fseq == 4 || fseq == -4) {
		// next a simple sort 
		for (int k = 0; k<_nDmd/2-1; k++) {
			for (int j=k+1;j<_nDmd/2; j++) {
				if ( (hopA[seq[k*2]]-hopA[seq[j*2]])*fseq <0){
					n1 = seq[k*2];
					seq[k*2] = seq[j*2];
					seq[j*2] = n1;
				}
			}
		}
	}
	//sort according to demand BW -- M(dmd)
	if (fseq == 2 || fseq == -2 || 
		fseq == 3 || fseq == -3 ||
		fseq == 4 || fseq == -4 ) 
		for (int k = 0; k<_nDmd/2-1; k++) {
			for (int j=k+1;j<_nDmd/2; j++) {
				if (fseq == 4 || fseq == -4 ) {
					if ( (M(seq[k*2])-M(seq[j*2]))*fseq > 0){
						n1 = seq[k*2];
						seq[k*2] = seq[j*2];
						seq[j*2] = n1;
					}
				} else if ( (M(seq[k*2])-M(seq[j*2]))*fseq <0){
					n1 = seq[k*2];
					seq[k*2] = seq[j*2];
					seq[j*2] = n1;
				}
			}
		}

	//	return 0;
	return fseq;
}

int ssr::global(double * totalCost)
{
	int k, m, dmd;
	int alldone;

	// decide the weight based on hop count.  Maybe used for seq selection
	int temp = 0;
	for (k = 0; k<_nDmd/2; k++) {
		seq[k*2] = k*2; seq[k*2+1]=k*2;
		Cflag[k*2] = 0; Cflag[k*2+1] = 0;
	}
	Ccurrent = 1;

	// decide sequence randomely using different seeds
	reshuffle();

	// resort the flows according to their working path hop count and nodes
	resort(fseq);  

	// next, do iterations several times for flows to update backup paths
	thisiter = 0;
	do {
		alldone = 0;
		reshuffle(); 
		resort(fseq); 
		for (k = 0; k<_nDmd/2; k++) {
			dmd = seq[k*2];
			// update backup path only when C changed recently.
			if ( (unsigned int)(Ccurrent - Cflag[dmd]) < (UINT_MAX-1)/2){
				if (pathr == 0 || pathr == 3 || pathr == 4)
					alldone += findbackup(dmd);
#ifdef FAILDENP
				else if (pathr ==1 || pathr ==2 )
					alldone += findbackupFD(dmd);
#endif
				else 
					assert(1);
			}
			else
				temp = 0;
		}

		thisiter ++;
		
	// finish when alldone and iteration count reached
	}while (thisiter < nIter && alldone > 0 ) ; 

#ifdef FAILDENP
	// for Spare reservation...
	spareBW = 0;
	spareLBW = 0;
	*totalCost = 0;

	if (cst != 3){  
		//getrsvBW();
		for (k=0; k<_nLink*2; k++){
			spareBW += (rsvBW[1][k]*unitcost[k]);
			*totalCost += NLC[(rsvBW[1][k]+rsvBW[0][k])*unitcost[k]] ;
		}
		if (multi==3){
			for (k=0; k<_nFail; k++){
				spareLBW += (rsvBW[2][k]);
			}
			//return spareLBW;
		}
	} else{// No share case
		for (k=0; k<_nDmd; k++)
			for (m=0; m<_nLink*2; m++){
				spareBW += (B[0][k][m]*M(k) );
				*totalCost += NLC[(B[0][k][m]+A[k][m])*M(k)] ;
			}
	}

	// for pathr == 3 again
	if (pathr == 3 || pathr == 4) {
		for (dmd = 0; dmd<_nDmd; dmd++) 
		  for (k=1;k<_nFail; k++) {
			memcpy(B[k][dmd],B[0][dmd],sizeof(char)*MAX_LINK*2);
			memcpy(route[k][dmd],route[0][dmd],sizeof(char)*MAX_NODE);
			succ2rt(dmd, rtnode[k][dmd], route[k][dmd]);
			hop[k][dmd] = hop[0][dmd];
		  }
	   //pathr = 1;
	   thisiter = 0;
	   do {
		alldone = 0;
		for (k = 0; k<_nDmd/2; k++) {
			dmd = seq[k*2];
			// update backup path only when C changed recently.
			if ( (unsigned int)(Ccurrent - Cflag[dmd]) < (UINT_MAX-1)/2){
					alldone += findbackupFD(dmd);
			}
			else
				temp = 0;
		}

		thisiter ++;
		
	// finish when alldone and iteration count reached
		}while (thisiter < nIter && alldone > 0 ) ; 
	}
#endif

	// for Spare reservation...
	spareBW = 0;
	spareLBW = 0;
	*totalCost = 0;
	if (cst != 3){  
		//getrsvBW();
		for (k=0; k<_nLink*2; k++){
			spareBW += (rsvBW[1][k]*unitcost[k]);
			*totalCost += NLC[(rsvBW[1][k]+rsvBW[0][k])*unitcost[k]] ;
		}
		if (multi==3){
			for (k=0; k<_nFail; k++){
				spareLBW += (rsvBW[2][k]);
			}
			return spareLBW;
		}
	} else{// No share case
		for (k=0; k<_nDmd; k++)
			for (m=0; m<_nLink*2; m++){
				spareBW += (B[0][k][m]*M(k) );
				*totalCost += NLC[(B[0][k][m]+A[k][m])*M(k)] ;
			}
	}

	//cout << thisiter << "\t";
	//cout << spareBW << endl;
	//cout.flush();

	return spareBW;
}

double ssr::comp_availability (void)
{
	double avg_avail = 0.0;

	int i;
	for (i=0; i<_nDmd; i++)
	{
		avg_avail += comp_avail_demand (i);
	}

	avg_avail /= _nDmd;

	return (avg_avail);

}

double ssr::comp_avail_demand(int i)
{
	double dmd_avail = 0.0;
	// for demand i, get the length of its backup route;
	// then get how long of its working route
	// store them, 
	// then compute how much degrade due to spare capacity sharing.
	dmd_avail = prob;

	return dmd_avail;
}